library(parallel)
library(doMC)
library(DHARMa)
library(broom.mixed)
library(MuMIn)
library(mgcv)
library(bbmle)
library(lme4)
library(merTools)
library(magrittr)
library(gridExtra)
library(ggplot2)
library(tidyr)
library(purrr)
library(plyr)
library(dplyr)1 - Variables
table 1: Variáveis utilizadas na descrição estatística
| code | description | class | range |
|---|---|---|---|
| n_nRef | number of SADs not refuted by the KS bootstrap test with critical p-value 0.05 | continuous | (0 ; 100) |
| diffS_mean | diffS = (S_MN - S_obs) / S_obs, diffS_mean = mean(diffS); S_MN = species richness estimated by MN | continuous | (-0.977 ; 4.78) |
| U_med | average of 10 replicates of the speciation rate estimated by MNEE | continuous | (8.860e-5 ; 4.221e-2) |
| p | proportion of tree cover | continuous | (0.013 ; 0.961) |
| MN | Neutral Model (EE = spatial explicit; EI = spatial implict) | category | 2 levels |
| d | mean dispersal distance (meters) | continuous | (0.456 ; 19.183) |
| k | proportion of propagules until the nearest neighborhood seq(0.99 : 0.05) | category | 20 levels |
| d_Lplot | mean dispersal distance / side of the sample area | continuous | (0.02 ; 0.192) |
| S_obs | observed species richness | integer | (26 ; 458) |
| Ntotal | number of individuals in the sample area | integer | (649 ; 12105) |
| SiteCode | sample area code | category | 103 levels |
Figure 1 Possible Predictor Variables
2 - Number of unrefuted SADs
Figura 2.1 número de SADs não refutadas ~ p * k * MN. A linha azul é uma estimativa baseada em ‘loess’.
Figura 2.2 número de SADs não refutadas ~ (d * MN|Site ~ p). Site está ordenado pelo valor de p.
2.1 GLMM binomial
Modelo cheio 1:
linear term
- ~ (p + p^2) * (d + d^2) * MN; if : var_dispersao = contiguous
- ~ ( p + p^2) * k * MN; if : var_dispersao = category
random term
- 1|SiteCode
- MN|SiteCode
- (var_dispersao + var_dispersao^2)*MN|SiteCode, if : var = dispersao : contiguous
# dados
df_md <- df_resultados %>%
select(n_nRef,p,SiteCode,MN,p.z,k,d.z,d_Lplot.z)
df_md[,6:8] <- apply(df_md[,6:8],2,as.character)
df_md %<>%
gather(key = var_dispersao, value = value_dispersao, k:d_Lplot.z)
# funcao ajuste dos modelos
f_md <- function(dados){
var_dispersao <- unique(dados$var_dispersao)
if(var_dispersao == "k"){
l_md <- vector("list",2)
names(l_md) <- c(paste0(var_dispersao," 1|Site"),
paste0(var_dispersao," MN|Site"))
dados$value_dispersao <- factor(dados$value_dispersao,levels = unique(dados$value_dispersao)[20:1])
l_md[[1]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * value_dispersao * MN +
(1|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
l_md[[2]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * value_dispersao * MN +
(MN|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
}else{
dados$value_dispersao <- as.numeric(dados$value_dispersao)
l_md <- vector("list",3)
names(l_md) <- c(paste0(var_dispersao," 1|Site"),
paste0(var_dispersao," MN|Site"),
paste0(var_dispersao," (",var_dispersao,"+",var_dispersao,"^2)*MN|Site"))
l_md[[1]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * (value_dispersao + I(value_dispersao^2)) * MN +
(1|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
l_md[[2]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * (value_dispersao + I(value_dispersao^2)) * MN +
(MN|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
l_md[[3]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * (value_dispersao + I(value_dispersao^2)) * MN +
( (value_dispersao + I(value_dispersao^2) ) * MN|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
}
return(l_md)
}
registerDoMC(3)
l_nRef__modeloCheio <- dlply(df_md,.(var_dispersao),f_md,.parallel = TRUE)
names(l_nRef__modeloCheio) <- NULL
l_nRef__modeloCheio <- do.call(c,l_nRef__modeloCheio)
#
f_warning <- function(glmm_object){
# glmm_object <- l_nRef__modeloCheio[[4]]
v_message <- glmm_object@optinfo$conv$lme4$messages %>%
as.character()
if(length(v_message)==0){v_message <- "OK"}
if(length(v_message)>1){v_message <- v_message[1] }
return(v_message)
}
df_auditoria.md <- ldply(l_nRef__modeloCheio,f_warning,.id="glmer") %>%
rename(warning_message = V1)Table 2.1 Modelos Cheios estimados e avisos de convergência
| glmer | warning_message |
|---|---|
| d_Lplot.z 1|Site | Model failed to converge with max|grad| = 0.00343548 (tol = 0.002, component 1) |
| d_Lplot.z MN|Site | Model failed to converge with max|grad| = 0.00836802 (tol = 0.002, component 1) |
| d_Lplot.z (d_Lplot.z+d_Lplot.z^2)*MN|Site | OK |
| d.z 1|Site | Model failed to converge with max|grad| = 40.1827 (tol = 0.002, component 1) |
| d.z MN|Site | Model failed to converge with max|grad| = 0.00746139 (tol = 0.002, component 1) |
| d.z (d.z+d.z^2)*MN|Site | OK |
| k 1|Site | OK |
| k MN|Site | OK |
allFit
i <- 1
names(l_nRef__modeloCheio[v_glmerUpdate])[i]
#
## update1:
md_allFit <- allFit(l_nRef__modeloCheio[v_glmerUpdate][[i]],maxfun = 1e5, parallel = 'multicore', ncpus = 3)1) d_Lplot 1|Site
7 optimizer(s) failed
2) d_Lplot MN|Site
7 optimizer(s) failed
3) d 1|Site
7 optimizer(s) failed
4) d MN|Site
7 optimizer(s) failed
Modelo cheio 2:
# dados
df_md <- df_resultados %>%
select(n_nRef,p,SiteCode,MN,p.z,k,d.z,d_Lplot.z)
df_md[,6:8] <- apply(df_md[,6:8],2,as.character)
df_md %<>%
gather(key = var_dispersao, value = value_dispersao, k:d_Lplot.z)
# funcao ajuste dos modelos
f_md <- function(dados){
var_dispersao <- unique(dados$var_dispersao)
if(var_dispersao == "k"){
l_md <- vector("list",2)
names(l_md) <- c(paste0(var_dispersao," 1|Site"),
paste0(var_dispersao," MN|Site"))
dados$value_dispersao <- factor(dados$value_dispersao,levels = unique(dados$value_dispersao)[20:1])
l_md[[1]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * value_dispersao * MN +
(1|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
l_md[[2]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * value_dispersao * MN +
(MN|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
}else{
dados$value_dispersao <- as.numeric(dados$value_dispersao)
l_md <- vector("list",4)
names(l_md) <- c(paste0(var_dispersao," 1|Site"),
paste0(var_dispersao," MN|Site"),
paste0(var_dispersao," ",var_dispersao,"*MN|Site"),
paste0(var_dispersao,"^2 (",var_dispersao,"+",var_dispersao,"^2)*MN|Site"))
l_md[[1]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * value_dispersao * MN +
(1|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
l_md[[2]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * value_dispersao * MN +
(MN|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
l_md[[3]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * value_dispersao * MN +
(value_dispersao * MN|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
l_md[[4]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * (value_dispersao + I(value_dispersao^2)) * MN +
( (value_dispersao + I(value_dispersao^2) ) * MN|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
}
return(l_md)
}
registerDoMC(3)
l_nRef__modeloCheio2 <- dlply(df_md,.(var_dispersao),f_md,.parallel = TRUE)
names(l_nRef__modeloCheio2) <- NULL
l_nRef__modeloCheio2 <- do.call(c,l_nRef__modeloCheio2)
#
f_warning <- function(glmm_object){
v_message <- glmm_object@optinfo$conv$lme4$messages %>%
as.character()
if(length(v_message)==0){v_message <- "OK"}
if(length(v_message)>1){v_message <- v_message[1] }
return(v_message)
}
df_auditoria.md <- ldply(l_nRef__modeloCheio2,f_warning,.id="glmer") %>%
rename(warning_message = V1)Table 2.2 Modelos Cheios estimados e avisos de convergência
| glmer | warning_message |
|---|---|
| d_Lplot.z 1|Site | Model failed to converge with max|grad| = 0.00284637 (tol = 0.002, component 1) |
| d_Lplot.z MN|Site | OK |
| d_Lplot.z d_Lplot.z*MN|Site | OK |
| d_Lplot.z^2 (d_Lplot.z+d_Lplot.z^2)*MN|Site | OK |
| d.z 1|Site | Model failed to converge with max|grad| = 0.00369235 (tol = 0.002, component 1) |
| d.z MN|Site | OK |
| d.z d.z*MN|Site | OK |
| d.z^2 (d.z+d.z^2)*MN|Site | OK |
| k 1|Site | OK |
| k MN|Site | OK |
allFit
d_Lplot.z 1|Site
7 optimizer(s) failed
d.z 1|Site
7 optimizer(s) failed
Comparação de Modelos Cheios:
Tabela 2.3 Comparação baseada em AICc dos modelos cheios
| GLMM | dAICc | df | weight |
|---|---|---|---|
| d^2 (d+d^2)*MN|Site | 0.0 | 39 | 1 |
| d_Lplot^2 (d_Lplot+d_Lplot^2)*MN|Site | 274.6 | 39 | <0.001 |
| d d*MN|Site | 38427.3 | 22 | <0.001 |
| d_Lplot d_Lplot*MN|Site | 38586.5 | 22 | <0.001 |
| k MN|Site | 74141.2 | 123 | <0.001 |
| d MN|Site | 95710.0 | 15 | <0.001 |
| d_Lplot MN|Site | 96279.6 | 15 | <0.001 |
| k 1|Site | 106158.0 | 121 | <0.001 |
Tabela 2.4 Coeficiente de Determinação Condicional e Marginal
## R2m R2c
## theoretical 0.3232656 0.9128366
## delta 0.3111874 0.8787302Diagnostico do modelo cheio mais plausível
Figura 2.3 Resíduos Quantílicos do modelo cheio plausível
Figura 2.4 Quantile-quantile plot random effects.
Figura 2.5 Predito e observado pelo modelo. Em vermelho linha y=x; em azul um modelo linear por sítio de amostragem.
Figura 2.6 Predito e observado para modelo cheio
Seleção da combinação de preditoras:
Modelo Global
md_nRef <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * (d.z + I(d.z^2)) * MN +
( (d.z + I(d.z^2) ) * MN|SiteCode),
family = "binomial",data=df_resultados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)),
na.action = "na.fail")Tabela 2.5 Tabela de Seleção de Modelos, deltaAICc<7
## Global model call: glmer(formula = cbind(n_nRef, 100 - n_nRef) ~ (p.z + I(p.z^2)) *
## (d.z + I(d.z^2)) * MN + ((d.z + I(d.z^2)) * MN | SiteCode),
## data = df_resultados, family = "binomial", control = glmerControl(optimizer = "bobyqa",
## optCtrl = list(maxfun = 2e+09)), na.action = "na.fail")
## ---
## Model selection table
## (Int) d.z d.z^2 MN p.z p.z^2 d.z:MN d.z:p.z d.z:I(p.z^2)
## 35712 2.296 -0.3911 0.03539 + 0.4903 -0.3671 + 0.5471
## 35840 2.382 -0.2154 0.03597 + 0.4569 -0.4532 + 0.4794 -0.17770
## 35696 1.929 -0.3912 0.03512 + 0.6350 + 0.5464
## 43904 2.296 -0.3940 0.03634 + 0.4741 -0.3668 + 0.6528
## 44032 2.383 -0.2168 0.03648 + 0.4409 -0.4545 + 0.5853 -0.17890
## 43888 1.929 -0.3939 0.03581 + 0.6170 + 0.6526
## 36736 2.265 -0.3906 0.07141 + 0.5017 -0.3356 + 0.5468
## 39808 2.435 -0.3912 0.03544 + 0.4367 -0.5064 + 0.5470
## 44928 2.264 -0.3932 0.07143 + 0.4864 -0.3345 + 0.6518
## 36864 2.363 -0.2319 0.04923 + 0.4622 -0.4336 + 0.4857 -0.16090
## 39936 2.411 -0.2225 0.03575 + 0.4459 -0.4824 + 0.4821 -0.17040
## 48000 2.436 -0.3939 0.03617 + 0.4201 -0.5077 + 0.6526
## 45056 2.364 -0.2318 0.04964 + 0.4481 -0.4352 + 0.5912 -0.16350
## 48128 2.414 -0.2237 0.03674 + 0.4291 -0.4855 + 0.5885 -0.17220
## 40832 2.380 -0.3905 0.06743 + 0.4577 -0.4514 + 0.5465
## 106368 2.501 -0.3906 0.08391 + 0.4095 -0.5720 + 0.5464
## 49024 2.381 -0.3932 0.06788 + 0.4410 -0.4525 + 0.6519
## 56320 2.389 -0.2016 0.03603 + 0.4525 -0.4603 + 0.4841 -0.19170
## 40960 2.392 -0.2376 0.04888 + 0.4503 -0.4610 + 0.4868 -0.15380
## 114560 2.502 -0.3933 0.08434 + 0.3934 -0.5732 + 0.6523
## 64512 2.423 -0.2323 0.03672 + 0.4255 -0.4949 + 0.5918 -0.16360
## 106496 2.478 -0.3213 0.07180 + 0.4188 -0.5488 + 0.5200 -0.07021
## 57344 2.384 -0.2154 0.04265 + 0.4555 -0.4550 + 0.4866 -0.17740
## 114688 2.479 -0.3211 0.07214 + 0.4020 -0.5501 + 0.6254 -0.07340
## 65536 2.408 -0.2846 0.06205 + 0.4309 -0.4792 + 0.6110 -0.11000
## I(d.z^2):MN I(d.z^2):p.z d.z^2):I(p.z^2 MN:p.z I(p.z^2):MN
## 35712 + -0.04297 +
## 35840 + -0.04304 +
## 35696 + -0.04152 +
## 43904 + -0.07235 +
## 44032 + -0.07353 +
## 43888 + -0.07158 +
## 36736 + -0.05651 -0.036240 +
## 39808 + -0.04324 + +
## 44928 + -0.08643 -0.035770 +
## 36864 + -0.04843 -0.013710 +
## 39936 + -0.04349 + +
## 48000 + -0.07313 + +
## 45056 + -0.07795 -0.013140 +
## 48128 + -0.07324 + +
## 40832 + -0.05550 -0.032450 + +
## 106368 + -0.06169 -0.048910 + +
## 49024 + -0.08515 -0.032170 + +
## 56320 + -0.04582 + +
## 40960 + -0.04801 -0.013600 + +
## 114560 + -0.09139 -0.048620 + +
## 64512 + -0.07326 + +
## 106496 + -0.05709 -0.036540 + +
## 57344 + -0.04806 -0.007001 + +
## 114688 + -0.08638 -0.035920 + +
## 65536 + -0.08312 -0.025920 + +
## d.z:MN:p.z d.z:I(p.z^2):MN I(d.z^2):MN:p.z I(d.z^2):I(p.z^2):MN df
## 35712 + 33
## 35840 + 34
## 35696 + 32
## 43904 + + 34
## 44032 + + 35
## 43888 + + 33
## 36736 + 34
## 39808 + 34
## 44928 + + 35
## 36864 + 35
## 39936 + 35
## 48000 + + 35
## 45056 + + 36
## 48128 + + 36
## 40832 + 35
## 106368 + + 36
## 49024 + + 36
## 56320 + + 36
## 40960 + 36
## 114560 + + + 37
## 64512 + + + 37
## 106496 + + 37
## 57344 + + 37
## 114688 + + + 38
## 65536 + + + 38
## logLik AICc delta weight
## 35712 -14179.02 28424.6 0.00 0.113
## 35840 -14178.06 28424.7 0.11 0.107
## 35696 -14180.31 28425.1 0.54 0.086
## 43904 -14178.29 28425.2 0.58 0.085
## 44032 -14177.32 28425.3 0.66 0.081
## 43888 -14179.58 28425.7 1.13 0.064
## 36736 -14178.68 28425.9 1.34 0.058
## 39808 -14178.82 28426.2 1.64 0.050
## 44928 -14177.96 28426.5 1.94 0.043
## 36864 -14178.03 28426.7 2.08 0.040
## 39936 -14178.05 28426.7 2.13 0.039
## 48000 -14178.09 28426.8 2.21 0.037
## 45056 -14177.28 28427.2 2.62 0.031
## 48128 -14177.31 28427.3 2.67 0.030
## 40832 -14178.55 28427.7 3.13 0.024
## 106368 -14177.79 28428.2 3.64 0.018
## 49024 -14177.83 28428.3 3.71 0.018
## 56320 -14178.00 28428.7 4.06 0.015
## 40960 -14178.04 28428.7 4.14 0.014
## 114560 -14177.06 28428.8 4.22 0.014
## 64512 -14177.30 28429.3 4.69 0.011
## 106496 -14177.72 28430.1 5.53 0.007
## 57344 -14177.99 28430.7 6.08 0.005
## 114688 -14177.01 28430.8 6.16 0.005
## 65536 -14177.20 28431.1 6.55 0.004
## Models ranked by AICc(x)
## Random terms (all models):
## '(d.z + I(d.z^2)) * MN | SiteCode'